{-# LANGUAGE ForeignFunctionInterface #-}

module Sim.Helper
  ( pair
  , unsafeSequenceIO
  , verifyPoint
  , verifyFoV
  , floor'
  , ergsToPhotons
  , quad
  , (<~)
  , parEval, mkChan
  , CacheSize, memo, memo2

  ) where

import   Control.Exception (evaluate)
import  Control.Concurrent 
import       Control.Monad ( replicateM )
import Control.Applicative ( (<$>) )
import           Data.List ( sortBy )
import       Data.Function ( on )
import          Data.Array ( Array, Ix, array, (!), bounds, range, inRange )
import    Data.VectorSpace ( (^/), (*^), (^-^) )
import    System.IO.Unsafe ( unsafeInterleaveIO )
import        NumInstances ()

type R = Double
type Z = Int
type CacheSize = ((Z,Z), (Z,Z))

foreign import ccall unsafe "math.h floor"
    c_floor :: R -> R

floor' :: R -> Z
{-# INLINE floor' #-}
floor' x = (truncate :: R -> Z) (c_floor x)

-- Втягивает функцию на множество пар.
pair :: (a -> b) -> (a, a) -> (b, b)
pair f = uncurry ( (,) `on` f )

--
unsafeSequenceIO :: [IO a] -> IO [a]
unsafeSequenceIO []     = return []
unsafeSequenceIO (x:xs) = do
   a <- x
   b <- unsafeInterleaveIO (unsafeSequenceIO xs)
   return (a:b)

-- Подтверждение или отрицание того, что существуют узел на указанной сетке,
-- координаты которого отличаются от координат указанной точки не более,
-- чем на 0.001.
verifyPoint :: (R,R) -> ((R,R), R) -> Bool
verifyPoint xy (zp,h) = a && b
   where 
     (a,b) = pair (<0.01) $ abs $ xy - uv

     uv = zp + (h *^ pair (fromIntegral.round) $ (xy - zp) ^/ h)
    

-- Подтверждение или отрицание того, что размер поля зрения меньше, чем
-- "размер" указанной сетки.
verifyFoV :: ((R,R), (R,R)) -> ((R,R), R) -> (Z,Z) -> Bool
verifyFoV fov (zp,h) (n,m) = not $ u < x || v < y || u' > x' || v' > y'
   where
     ((u,v), (u',v')) = fov
     ((x,y), (x',y')) = (zp, zp + h *^ (pair fromIntegral (n,m) - 1))

-- Пересчет эрги в фотоны
ergsToPhotons f = \x -> (f x) * x / hc
   where
    -- Произведение постоянной Планка на скорость  света [erg * А / sec^-1]
    hc = 6.62606896E-27 * 2.99792458E+18

-- Интегрирование функции методом трапеций.
quad :: (Fractional a, Enum a) => (a,a) -> a -> (a -> a) -> a
quad (a,b) d f = d * (c + e)
   where
      c = 0.5 * (f a + f b)
      e = sum $ [f x | x <- [a+d,a+2*d .. b-d]]

-- Трансформация массива в чистую функцию так, что вне диапазона задания
-- массива функция имеет значения равные граничным значениям этого массива.
(<~) :: (Ix i, Ix j) => Array (i,j) e -> (i,j) -> e
(<~) arr (i,j) = arr ! (index bi i, index bj j)
   where
      index (l,r) i | i < l = l
                    | i > r = r
                    | otherwise = i

      bi = (fst l, fst r)
      bj = (snd l, snd r)

      (l, r) = bounds arr

-- Мемозация (кэширование) функции
memo :: CacheSize -> ((Z,Z) -> b) -> (Z,Z) -> b
memo r f = memo'
  where
    memo' xy | r `inRange` xy = arr ! xy
             | otherwise      = f xy

    arr = array r [ (i, f i) | i <- range r ]

memo2 :: CacheSize -> ((Z,Z) -> (Z,Z) -> b) -> (Z,Z) -> (Z,Z) -> b
memo2 cb f =  memo cb (memo cb . f)


--


parEval :: (Ord k, Enum k, Num k) =>  [v] -> (Chan (k, v), Chan (k, v)) -> IO [v]
parEval [] _         = return []
parEval xs (inp,out) = do

  writeList2Chan inp $ zip [1..] xs

  es <- map snd . sortBy sortFunc <$> replicateM n (readChan out)

  return es

  where
    n = length xs

    sortFunc (a,_) (b,_) | a == b = EQ
                         | a >  b = GT
                         | a <  b = LT

mkChan :: Int -> IO (Chan (k, v), Chan (k, v))
mkChan n = do

--  n   <- getNumCapabilities

  inp <- newChan
  out <- newChan

  mapM_ (forkIO . evalFunc inp) $ replicate n out

  return (inp,out)

  where
    evalFunc inp out = do
       
        (i,x) <- readChan inp

        e <- evaluate x

        writeChan out (i,e)
        
        evalFunc inp out